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ABSTRACT 

Weak gravitational lensing provides a unique method to directly measure the distribu- 
tion of mass in the universe. Because the distortions induced by lensing in the shape of 
background galaxies are small, the measurement of weak lensing requires high preci- 
sion. Here, we present a new method for obtaining reliable weak shear measurements. It 
is based on the Shapelet basis function formalism of Refregier (2001), in which galaxy 
images are decomposed into several shape components, each providing independent 
estimates of the local shear. The formalism affords an efficient modelling and decon- 
volution of the Point Spread Function. Using the remarkable properties of Shapelets 
under distortions, we construct a simple, minimum variance estimator for the shear. 
We describe how we implement the method in practice, and test the method using 
realistic simulated images. We find our method to be stable and reliable for conditions 
analogous to ground-based surveys. Compared to earlier methods, our method has 
the advantages of being accurate, linear, mathematically well-defined, and optimally 
sensitive, since it uses the full shape information available for each galaxy. 

Key words: methods: data analysis; techniques: image processing; cosmology: ob- 
servations; dark matter; gravitational lensing; large-scale structure of Universe 



1 INTRODUCTION 

Weak gravitational lensing is a powerful method to map the 
mass of clusters of galaxies (for reviews see Mellier 1999; 
Bartelmann & Schneider 2000). Recently this technique has 
been extended to large scale structures by several groups 
(Wittman et al 2000; van Waerbeke et al 2000; Bacon, Re- 
fregier & Ellis 2000; Kaiser et al 2000; Maoli et al 2001; 
Rhodes, Refregier & Groth 2001; van Waerbeke et al 2001), 
and thus offers bright prospects for cosmology. 

Because the lensing effect is only of a few percent on 
large scales, a precise method for measuring the shear is re- 
quired. The widely used method of Kaiser, Squires & Broad- 
hurst (KSB, 1995) has several shortcomings in the context 
of upcoming weak lensing surveys (see also the early method 
of Bonnet & Mellier 1995). Firstly, it is not sufficiently ac- 
curate or unbiased to measure shears of a fraction of 1% (cf 
Bacon et al 2001a, Erben et al 2001). Secondly, KSB suf- 
fers from being mathematically ill-defined (cf Kaiser 2000, 
Kuijken 1999) for space- and ground-based PSFs. 

Thus, several new methods have been proposed for 
measuring weak lensing. Rhodes, Refregier & Groth (2000) 
proposed a variant of KSB tuned to the analysis of Hub- 
ble Space Telescope (HST) images. Kaiser (2000) advo- 
cates a new shear estimator based on second moments of 
galaxy shapes, which correctly deals with desmearing. Kui- 
jken (1999) presents a method which fits several gaussian 
profiles in order to account for both smearing and shearing. 



Here, we present an independent approach, based on 
formalism introduced in Refregier (2001, Paper I). In this 
approach, galaxy images are linearly decomposed into a 
series of orthogonal basis functions of different shapes, or 
'shapelets'. Because of the remarkable properties of the ba- 
sis functions, they are particularly well suited for shear es- 
timation. We describe a process for deconvolving the Point 
Spread Function (PSF), and for estimating shear in a well- 
defined manner using the shapelets. We show how our 
method can be implemented in practice and test its reli- 
ability using realistic numerical simulations. Compared to 
earlier methods, our method has the advantages of being 
accurate, linear, mathematically well-defined, and optimally 
sensitive, since it uses the full shape information available 
for each galaxy. This will be particularly important to fully 
exploit future space-based weak lensing surveys with HST 
and the future SNAP mission (Perlmutter et al. 2001), and 
ground-based wide-field surveys such as those with Mega- 
cam (Bernardeau et al. 1997), VISTA (Taylor et al. 2001), 
DMT (Tyson et al. 2000) and WFHRI (Kaiser et al. 2000). 
An application of our method to interferometric images will 
be presented in Chang & Refregier (2001). 

The paper is organised as follows. In we collect to- 
gether the necessary formalism for our shapelet description 
of galaxies. In §^ we describe the shapelet (de) convolution 
matrix and explain how to correct object shapes for the 
PSF. In we present the shear estimators which exist for 
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each shapelet coefficient for a galaxy, and we describe how 
to combine these estimators to maximise the signal. In §H 
we turn to the practical implementation of the method, and 
in we describe the results of testing the method on sim- 
ulated data. Our conclusions are summarised in 



2 SHAPELET FORMALISM 

We begin by summarising the necessary formalism from Pa- 
per I for a description of galaxies in our basis set. A galaxy 
with intensity /(x), can be decomposed into our basis func- 
tions -Bn(x; P) as 



/(x)=^/„Bn(x;/3). 



(1) 



where x — {xi,X2) and n = (ni,n2). The 2-dimensional 
cartesian basis functions can be written as Bn{yi;P) = 
-B„i (xi; /3)_Bn2 (a;2; in terms of the 1-dimensional basis 
functions 



Bn{x-P) 



2"7r^n!/3 



e 2^ 



(2) 



where H„{x) is a Hermite polynomial of order n. The pa- 
rameter /3 is a characteristic scale, which is typically chosen 
to be close to the radius of the object (see discussion in 
Paper I, § 3.1, and ^ below). 

Because these basis functions, or 'shapelets' form a com- 
plete orthonormal set, the coefficients /„ can be found using 



d^x /(x)Bn(x;/3). 



(3) 



This decomposition provides an excellent and efficient de- 
scription of galaxy images in practice (see Paper I for exam- 
ples of decomposition and recomposition of galaxy images) . 
Note that the ni + n2 = 2 shapelet coefficients are exactly 
equal to the gaussian-weighted quadrupole moments used 
in the KSB method. Our shapelet method thus, in a sense, 
generalises the KSB method to use all available multipole 
moments. 



3 DECONVOLUTION OF THE 
POINT-SPREAD FUNCTION 

Our first concern in providing a measure of the shear is 
to remove the the effect of the PSF convolution or 'smear- 
ing', which acts upon galaxy images. The PSF is generally 
anisotropic and results from atmosphere turbulence or 'see- 
ing' (for ground-based observations) , tracking errors, imper- 
fect optics, etc. Since typical PSF ellipticities are of order 
10% while the sought-for shear signal is of order 1%, an 
accurate correction for the PSF is vital. Here, we describe 
the convolution formalism for shapelets in detail, and then 
discuss how it can be used to deconvolve the PSF. 



3.1 Convolution Formalism 

We now show how our basis functions behave under convolu- 
tions. For this purpose, let us consider a galaxy of intensity 
/(x) observed with an instrument with a PSF (^(x). The 
observed image /i(x) is given by the convolution 



Table 1. First few componentst of the normalised 3-product 
integral L;„„(a, fe, c) 



iooo = 1 

i002 = -2 + 2c2 

Loii = 2c6 

Lo22 = 4 - 4fe2 - 4c2 -I- 12fe2c2 

Lii2 = -4afe -I- Uabc^ 

Loi3 = -I2bc+12bc^ 

Lo04 = 12 - 24c2 + 12c'' 

Loo6 = -120 + 360c2 - seOc-i + 120c^ 

t Other components can be obtained by symmetry (eg. 
^020 = ~2 + 2b^); Components with odd I + m + n vanish. 



/j(x) = if * 3)(x) = / d'x' /(x - x')ff(x') 



(4) 



Using Equation (|3|), we can first decompose each of these 
three functions into their shapelet coefficients /n, gn, hn 
with shapelet scales a, j3 and 7, respectively. As discussed 
in Paper I, the convolved coefficients are related to the un- 
convolved ones by 



/in 



E 

ml 



(5) 



where Cnmi(7, ct, P) is the 2-dimensional convolution tensor, 
which can be written in terms of the 1-dimensional convo- 
lution tensor C„mi(7, a, /3) as 



Cnml(7, a, /3) = Cnimiii (7, «, P)Cn2m2l-2 {l , «, P) ■ 



(6) 



Using the invariance of the basis functions under Fourier 
transform (see Paper I, § 2.2), we find that this latter tensor 
is given by 

C„,„,(7,a,/3) = (2^)^(-l)"i"+™+'Biti,(7-\«-\/3-'),(7) 
where -^^'^^(ai, a2, as) is defined as 

Bl^„{ai,a2,a:i) = / dx Bi{x,ai)B^{x,a2)B„{x,a3). (8) 



B 



(3) 



We now seek to evaluate the key 3-product integral 
For this purpose, we first rewrite it as 



2 TT 2 m\n\l\aia2a3 

V ai a2 a-i 



(9) 



where u = a-^ + "2 + "^3 ^^'^ where we have defined 



Limn{a,b,c) = —= 



dx e"" Hi{ax)Hm{bx)H„{cx). (10) 



By parity this integral vanishes ii m + n + l is odd. By using 
the relation between H„-i{x) and H'nix) and by integrating 
by parts, one can show that this integral obeys the recur- 
rence relation 

Li+\^m,n ~ 2/(a^ — l)L;_i_m_„-|-2ma6Li_m_i,„-(-2nacL;_,„_„_i ,( 

and similarly for L;,m+i,n and Li,m,n+i. This and the fact 
that Looo = 1 can be used to conveniently evaluate Bf ^ 
The first few components of Limn{a,h,c) are listed in Ta- 
ble §. 
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We therefore have derived a fuUy analytical expression 
for the convolution of two functions within the shapelet for- 
malism. Indeed, to convolve an object / with an object g, 
one can decompose / and g into shapelets, calculate the 
C tensor using Equations (p|)-(p^, and then apply Equa- 
tion Q to obtain the convolved object h. Note that the 
above recurrence relation (Eq. is also very useful for 

other problems, such as deprojection and Poisson noise, in 
which the 3-product integral B^^^ naturally arises (see Pa- 
per I). 



3.2 Deconvolution Method 

We are now in a position to develop a deconvolution method. 
We are seeking an approach that allows us to decompose the 
PSF (measured from the stars) and the smeared galaxy in 
question, and obtain directly the deconvolved coefficients 
which can then be the input to our shear estimators. Be- 
cause of the above properties of the basis functions under 
convolution, this is a simple matter of linear algebra. 

Let us assume that the PSF (?(x) has been measured 
(typically from stellar images) and decomposed into its 
shapelet coefficients as described above. It is then con- 
venient to combine these coefficients with the convolution 
tensor in Equation (Bl) and write 




E 



Pnmfn 



(12) 



where we have defined Pnm = CnmiSi- By arranging the 
coefficients n = (711,712) into a vector, Pnm can be consid- 
ered as a matrix, which we call the 'PSF matrix'. 

Our goal is to recover the unconvolved coefficients /„ 
from the observed (convolved) coefficients hn. One way to 
achieve this is to attempt to invert the PSF matrix. In Pa- 
per I, however, it was shown that convolution amounts to 
a projection of the high-order shapelet states onto states of 
lower order (see also ilUustration in Figures |^ and ^ . This 
is expected, since the high-order modes have high frequency 
oscillations and are thus smeared out by the convolution. 
As a result, the PSF matrix has typically small high-order 
entries and is thus not invertible as is. On the other hand, if 
we restrict the PSF matrix to entries of sufficiently low or- 
ders, it will indeed be invertible. This amounts to giving up 
on the recovery of high order information, which has been 
destroyed by convolution. In we discuss how we choose 
a, 13 and 7, and the maximum order of recovery for the con- 
volution in practice. 

After this restriction to low order, we can thus invert 
the PSF matrix and obtain 



/m — Pmn^n- 



(13) 



This provides an estimate for the (low order) coefficients of 
the unsmeared object /. 

As Kuijken (2000) suggested, another approach consists 
in fitting the observed galaxy coefficients h = {/in} for the 
deconvolved coefficients f = {./n}. This can be done by min- 
imizing 



X 



(h-Pf)^V~^(h-Pf) 



(14) 



with respect to the model parameters f given the data vector 
h. Here, Kim = cov(/in, ftm) is the covariance of the observed 




1 00 1 50 



Figure 1. Example of deconvolution by shapelet decomposi- 
tion: (a) Image before smearing, (b) Convolving kernel (PSF), 
(c) Smeared image, (d) Recovered image from the shapelet ma- 
trix inversion. 
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Figure 2. Coefficient matrices at different stages of tlie deconvo- 
lution sfiown in Figure 1: (a) original image, (b) convolved image, 
(c) deconvolved image. 



coefficients hn resulting from noise in the observed image. 
As discussed in Paper I, this can easily be evaluated from the 
properties of the pixel noise. For instance, it is proportional 
to the identity matrix in the case of white noise. Since this 
model is linear in the model parameters f , the best fit has 
the simple analytic solution (see eg. Lupton 1993) 



(15) 



with a covariance error matrix Wnm = cov(/n, /m) given by 
W = (P"^V~'^P)~^. These analytic solutions make the best 
fit parameters and errors efficient to evaluate in practice. 

The latter scheme is potentially more robust numeri- 
cally since it does not require the direct inversion of the 
PSF matrix. On the other hand, it is more complicated and 
requires knowledge of the noise properties of the observed 
image. In particular, the procedure is, strictly speaking, 
only valid in the case of gaussian noise, a condition never 
rigorously met in optical CCD images. We have tried both 
deconvolution schemes on various object and PSF shapes 
and have found that they both perform equally well (as 
long as the PSF matrix is restricted to low-order as dis- 
cussed above). For the remainder of this paper, we will use 
the direct inversion scheme, i.e. Equation (|l3|). 

Figure |l| gives an illustration of the procedure. A galaxy 
(panel 1) is smeared by a complicated PSF (panel 2). The 
resulting object is shown in panel 3. The deconvolved galaxy 
image obtained using our direct inversion scheme (using 
?imax = 9, /3 = 15 and a — 'y — 20 pixels) is shown on 
panel 4. It is very close to panel 1, showing that our method 
is successful in recovering the original image. A deconvolu- 
tion using the method yields very similar results. Figure 
shows the shapelet coefficient matrix at different stages 
of the procedure. Notice the projection of high-order coeffi- 
cients for the galaxy (panel 1) into lower-order coefficients 
after convolution (panel 2) . The recovered coefficients (panel 
3) are very close to the original coefficients (panel 1). 



4 MEASURE OF THE SHEAR 

Now that we have a method to correct for the PSF, we turn 
to the problem of measuring the shear from an ensemble of 
galaxies. We start by presenting the formalism used to de- 
scribe the action of shear with shapelets. We then construct 
shear estimators for each shapelet coefficient and combine 
them to derive a minimum variance estimator. 



4.1 Shear Matrix 

Let us consider a galaxy with an unlensed intensity ,/"(x). 
In our formalism (see Paper I), the lensed intensity after a 
weak shear ji is written as 



I 



f'~{i + jA)f, 



(16) 



to first order in the shear, where Si is the shear operator. 
If we decompose these intensities into our basis functions 
i3n(x, /3) (Eq. Q), this can be expressed as a relation be- 
tween the lensed and the unlensed coefficients: 



(17) 
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where Simn = / d^a;i3n(x)S'ii3m(x) is the shear matrix. As 
presented in Paper I, our basis functions are also the eigen- 
functions for the Quantum Harmonic Oscillator (QHO), 
thus allowing us to use the powerful formalism developed 
for this problem. In particular, the shear operators can be 
written as 



xlO" 



Si 



At2 



-2 



§2 = a,\a^ 



— CL\a2, (18) 

where a\ and di are the raising and lowering operators, for 
each dimension i = 1,2. They operate as 



&lBn 



niBn 



= Vril + Ii3ni + l,n2,(19) 

and similarly for a2 and a\. The shear matrices are thus 
simple to evaluate in this way, and are very sparse since 
they involve the mixing of only a few modes. 

An illustration of the action of the shear matrix on the 
first few shapelet states can be found in Paper I. A more re- 
alistic example is shown on Figure ^. Here, we have used the 
shear matrix to shear a galaxy (panel 1) found in the Hub- 
ble Deep Field (Williams et al. 1996) by 20%. The resulting 
sheared image (panel 2) is virtually indistinguishable from 
the same image sheared directly in real space (panel 3). The 
action of shear on the shapelet coefficients is illustrated in 
figure ^ The average shapelet coefficients for a population 
of galaxies in one of our simulations (see ^ is shown in the 
upper panel. The difference between the coefficients before 
and after shears of 20% in the 71 and 72 directions are shown 
in the middle and lower panels, respectively. Note the change 
in even-even coefficient components associated with a shear 
in the 71 direction, and the change in odd-odd coefficients 
for a shear in the 72 direction (see below for a discussion of 
this effect). 

4.2 Shear Estimator 

Our goal is to find an estimator 7^ for the shear 7^. We re- 
quire that this estimator be unbiased, i.e. that (7^) = ji, 
when averaged over a population of galaxies which are ran- 
domly oriented before lensing. To limit the impact of noise 
on our estimator, we also require it to be linear in the 
sheared coefficients f^. 

To construct such an estimator, we first notice that the 
average shapelet coefficients (/n), before lensing, must be 
rotationally invariant. This must be so since the ensemble 
of randomly-oriented galaxies does not have a preferred di- 
rection. Using the polar shapelets discussed in §5 of Paper 
I, it is easy to show that this will be satisfied only if (/n) 
vanishes when ni and/or n2 is odd. This can be verified by 
inspecting the top panel of Figure ^ which shows that the 
only non-zero unlensed coefficients are the even-even ones, 
in a simulated galaxy ensemble (see 

The coefficients after lensing will no longer have this 
symmetry, since the shear introduces a preferred direction. 
In our formalism, this results from the mixing induced by 
the shear operator (see Eq. ^|). From Equation (p^, it 
is easy to see that a 71-shear mixes even-even states along 
the vertical and horizontal axes in the ni — 712 shapelet co- 
efficient plane. On the other hand, a 72-shear mixes states 
diagonally on this plane. As a result, 71 only affects states 
with Til and 722 even, while 72 only those with ni and n2 

© 2001 RAS, MNRAS 000, Flfl 
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Figure 3. Example of the action of the shear matrix on a galaxy 
found in the Hubble Deep Field. The original image (panel 1) 
is sheared by 20% by the shear matrix (panel 2). The resulting 
image is almost indistinguishable from that of the same galaxy 
sheared directly in real space (panel 3). 
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Figure 4. Action of sliear on the average sliapelet coefficients 
for a galaxy population in our simulation: (a) Initial simulation, 
(b) difference between initial and sheared simulation coefficients 
= n.9V (r^ .sifl.me fl-ci ^^h^ for = n.9. 



odd. (States with ni even and 712 odd, or vice- versa, are left 
unchanged by shear). This convenient fact means that the 
two shear components are decoupled in shapelet space (see 
figure ^). We can therefore construct independent estima- 
tors for each component. It is easy to show, that the only 
estimators satisfying the above conditions are 



7ln 



72n = 



5'lnm {fm} 



ni , n2 even 



ni, n2 odd. 



(20) 



where only even-even (odd-odd) coefflcients are used for 7in 
(72n). As before, the brackets denote an average over the 
galaxy ensemble and can be estimated using a large region 
where the mean shear is effectively 0. This provides us with 
a shear estimator 7^1, for every (appropriate) shapelet coef- 
ficient of each galaxy. 

We now seek to combine these estimators in a manner 
which maximises the shear signal. For this purpose, we can 
construct a combined shear estimator of the form 



7i 



En 



(21) 



where the weights win {w2n) are set to zero when n is not 
even-even (odd-odd). By construction, 7^ is still linear and, 
thanks to the denominator, guaranteed to be unbiased. The 
weights then need to be chosen to minimise the uncertainty 
in 7i . To find the optimal weights, we consider the covariance 
matrix between the individual estimators 



COv(7in, 7im), 



(22) 



which can be computed by averaging over the (unsheared) 
galaxy ensemble. It is easy to show that the variance var(7i) 
of the combined estimator is minimized when 



(23) 



With this optimal choice for weights uiin, the estimator vari- 
ance reduces to 



r(70 = I X! ^'"'^ 

\ nm 



(24) 



We have thus derived an optimal shear estimator 
(Eq. |2^) for each galaxy. It can then be averaged over all 
galaxies in a region to provide a local estimate of the shear. 
The error in the resulting shear measurement can then be 
computed using Equation (|2^), or directly from the variance 
of the actually measured shear estimators. We now discuss 
how to implement the shear estimation in practice. 



5 IMPLEMENTATION 

The above formalism is fairly straightforward, and is easily 
implemented. Here we describe some important practical de- 
tails which are required for measuring shear from real data. 

The objects in an image are first detected and cat- 
alogued using the publicly available software SExtractor 
(Bertin & Arnouts 1996). Each object is then cut out from 
the image, choosing a square centred upon the SExtractor 
centroid, extending by 5 times the SExtractor semi-major 
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axis found for the object. This ensures that the box size 
is much larger than the object, as is necessary for the or- 
thogonality of the basis functions. The median background 
around each of these objects is then subtracted. 

Stars are selected, either from a magnitude-radius plot 
or from the SExtractor neural network classifier. The stars 
are decomposed into shapelet coefficients, choosing /3 to be 
0.7 times the semi-major axis found by SExtractor (this 
choice is convenient as it leads to a compact PSF repre- 
sentation in shapelet space) . We iterate to find the best cen- 
troid for decomposition by decomposing a first time, finding 
a new centroid using Equation (27) in Paper I, then decom- 
posing a second time about the new centroid. For this paper, 
we set the maximum order of the stellar decomposition to 
?^max = 6. Each stellar shapelet component is divided by the 
stellar flux calculated from Equation (26) in Paper I, thus 
yielding normalised stellar coefficients. 

We can then interpolate each stellar normalised coeffi- 
cient across the field using a 2-dimensional polynomial. This 
affords us with a model of the PSF at each point on the im- 
age, which we can then use to deconvolve the galaxies. In the 
simulation described below (§^, the PSF is constant across 
the field. In this paper, we thus simply average the coeffi- 
cients from all stars, and use this average coefficient set for 
desmearing at all positions in the field. 

We now decompose each galaxy according to Equa- 
tion (^), using rimax = 6 as before. We choose the shapelet 
scale for the galaxies to be fixed to the median SExtractor 
FWHM for the set of galaxies in question (see discussion 
of binning below). In the case of galaxies where this scale 
is close to the decomposition scale for stars (/3), we set the 
decomposition scale to 1.5/3. As for the stars, we iterate to 
obtain the best centroid about which to decompose. 

We then deconvolve each galaxy with the PSF model 
at that position, using Equation ([L3|). To make the matrix 
inversion stable (see we only keep a finite number of 
elements of the PSF matrix, here n^ax = 6. (We could al- 
ternatively apply the modelling of Eq. Ilq]). The best 
results are obtained by choosing the deconvolved shapelet 
scale a to be equal to 7. Indeed, a choice of a comparable 
to 7 uses all the available information, and does not prohibit 
us from recovering small objects. 

The next step is to bin the galaxies by magnitude and 
radius. This is to reduce the dispersion in the shapelet coeffi- 
cients, and therefore the overall uncertainty of the combined 
shear estimators. The choice of bin size depends on the size 
of the data set, as each bin should contain enough objects 
to afford reasonable signal-to-noise in the bin. For each bin 
we obtain shear estimators 7in from all (appropriate) coeffi- 
cients n, using Equation (|20[). We then calculate the weight 
Win for each shear estimator according to Equation (p3|) , and 
compute the combined estimator for this galaxy (Eq. pif). 

Next we perform an (unweighted) average of the shear 
estimators for all galaxies within a given magnitude-radius 
bin. We can finally perform a weighted average of the shear 
estimators from all bins, using the variance within each bin 
as the inverse weight; we can even weight with magnitude 
in order to optimise the weak-lensing signal from a given 
redshift interval. 



6 SIMULATIONS 

We now describe our initial tests of this new method using 
simulations. Our goal is to verify that our method can re- 
cover shear which we impose upon simulated images with 
realistic observational properties. A further paper (Bacon 
et al 2001b) will describe a range of ground- and space- 
based applications, and compare results obtained with KSB 
and shapelets for simulated and real data. Here we restrict 
ourselves to ground-based applications, while recognising 
that the preservation of more than second moments with 
shapelets will have greatest initial impact on shear mea- 
surements from space. 

The simulations are based on those described in Bacon 
et al. (2001a). Full details are contained in that paper; here 
we briefly summarise the relevant features of the simula- 
tions. We create realistic simulations of flelds of galaxies ob- 
served by a typical ground-based 4m telescope (the William 
Herschel Telescope), with appropriate magnitudes, counts, 
diameters and ellipticities for stars and galaxies. We include 
an appropriate range of seeing, input shear, and tracking 
errors, and set the simulations in the context of the appro- 
priate pixel scale. 

We obtain the required galaxy statistics via the the re- 
solved (0.1 arcsec seeing) image statistics of the Groth Strip 
(Groth et al. 1994), a deep (7 ~ 26) survey observed by 
the Hubble Space Telescope (cf Bacon et al 2001a). Since 
the HST PSF is much smaller than that typical for WHT 
(0.7"), the Groth Strip provides effectively unsmeared el- 
lipticities and diameters. The Groth Strip contains approx. 
10,000 galaxies in a 108 arcmin^ area. We utilise Ebbels' 
(1998) SExtractor catalogue obtained from the strip, which 
contains magnitude, diameter, and ellipticity for each ob- 
ject. 

We flt a model to the multi-dimensional probability dis- 
tribution of galaxy properties (eg morphology, ellipticity, 
diameter, magnitude) in this catalogue. With the result- 
ing model, we can draw a statistically similar catalogue of 
galaxies with a realistic distribution of these properties, via 
Monte Carlo selection. We spatially distribute the objects 
randomly over a CCD frame, and add stars of appropriate 
magnitude distribution modelled from WHT data. 

It is now a simple matter to shear the galaxies in our 
catalogue. We flrst calculate the change in each object's el- 
lipticity ei which lensing will induce. To first order in the 
shear, the ellipticity transforms as (see Rhodes et al 2000) 



e- = e, 2{S,j + eiej)7j 



(25) 



Similarly, the magnitude of the sheared object is related to 
that of the initial object by 

m' = m + 2.51og(l-7^) (26) 

Finally, we compute the sheared semi-major axis size. In 
order to do this, we define R, A and J matrices by 



R = 



- sin (/) cos 4> 







(27) 



(28) 



J = R^AR (29) 
where (j) is the position angle of the object in question, a is 
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its semi-major axis, and b its semi-minor axis. Then, if we 
define a shear matrix $ as 

<i>=("" ) (30) 

we can find the new semi-major axis a' using 

Jk.l = Jk.l + JkA^lA + Jl,i^k.i (31) 

a" = 0.5 (^J^,o + J(,i + ^j{Jio-JliY+^Ji^^ ■ (32) 

These relations provide us with the sheared object cata- 
logue. 

For the purposes of this paper, we ran a set of simula- 
tions with shears in both 71 and 72 directions, ranging from 
zero to 5%. This affords a check of the shapelet method in 
the weak shear regime. We set a shear which is uniform over 
a given field. This is adequate for our initial tests, in which 
we are interested in recovering a mean shear across a whole 
field. 

We model tracking errors by giving an anisotropic PSF 
to the fields. Stellar ellipticities are chosen as uniform across 
a given field, with the variance of the ellipticity from field 
to field equal to (Te=0.05. 

We produce simulated images from the catalogues with 
the IRAF artdata package. This draws stars and galaxies 
from the catalogues with specified magnitude, diameter, el- 
lipticity, morphology (de Vaucouleurs or exponential) and 
position. Telescope-specific details are included: telescope 
throughput, anisotropic PSF (with Moffatt profile, and see- 
ing chosen to be 0.6"), pixellisation (0.24" per pixel), Poisson 
and read noise, sky background and gain are all realistically 
modelled. Examples of the resulting images can be found in 
Bacon et al (2001a). 

After image realisation, we apply our shear measure- 
ment algorithm to the simulated images. We train the esti- 
mators (i.e. obtain the necessary {/„) for Eq. [^) with an 
unsheared set of galaxies with the magnitude-radius distri- 
bution and intrinsic properties appropriate for our selected 
cell. For our testing purposes, we examine the shear in one 
narrow magnitude-radius bin only, centred on m = 22.5 and 
a = 2.0 pixels. 

We compare the input shear for our simulated fields 
with the shear estimates obtained by the shapelets method. 
The results are shown in figure I for the 0-5% shear sim- 
ulations described above. Notice that the output shear is 
linearly related to input shear, with a slope consistent with 
1. With linear regression we indeed find 7^"* = 0.977i°; and 
72^* = l-0072"i, with standard errors on the slope of 0.04 in 
each case. 

We conducted a further set of simulations to exam- 
ine the effect of seeing (PSF size) on our recovery. A set 
of simulations with identical shear (71,72) = (0.00,0.05) 
were produced, with seeing values ranging from 0.2" to 0.8". 
The recovery and noise are demonstrated in Figure ^. Our 
method clearly remains unbiased at all seeing values consid- 
ered. Since we examine the shear in a single, small magni- 
tude bin, the number density does not increase to improve 
the signal with decreased PSF, as is the case when all mag- 
nitudes are considered (see Bacon et al 2001). Nevertheless, 
a small decrease in noise is seen at small seeing values, show- 
ing that more information is recovered in this case. 
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-0.10 1/ I 

-0.10 -0.05 0.00 0.05 0.10 

7i,i„ 




Figure 5. Input shear vs recovered shear for a set of 11 simula- 
tions. 

These initial simulations demonstrate the effectiveness 
and utility of this shear measurement method. We have 
made additional tests to check the reliability of the method 
upon varying galaxy type, size and magnitude; the shear 
is recovered accurately in each case. Detailed discussion of 
these properties of the method on real and simulated data 
will be carried out in a further paper (Bacon et al 2001b). 



7 CONCLUSIONS 

In this paper, we have presented a new method for preci- 
sion measurements of weak lensing. After summarising the 
necessary formalism from Paper I, we described the means 
of convolving and deconvolving objects using shapelets. In 
shapelet space, this is a simple matrix operation, which can 
be used to correct for the PSF. We then presented shear 
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Figure 6. Results for set of simulations showing shear recovery 
with varying seeing. Note the excellent recovery at all seeing val- 
ues, with slight increase in noise with increasing PSF. 



estimators for shapelets. These are obtained from another 
simple matrix formalism, and are constructed to be unbi- 
ased and linear. Each shapelet coefficient (even-even and 
odd-odd) provides a shear estimator, thus using all the avail- 
able shape information for each object. We combined these 
individual estimators to construct a minimum variance esti- 
mator for the shear, ensuring that the shear signal is max- 
imised. 

The reliability of our shapelet method was then tested 
using simulations of realistic ground-based images. We find 
the method to be accurate and stable against variation of 
galaxy type, size of object, noise level, and PSF characteris- 
tics. In a future paper (Bacon et al 2001b), we will present 
detailed tests of our method based on extensive ground- and 
space- based image simulations. In particular, we will com- 
pare in detail the performance of our method against that 
of the commonly used KSB method. 

Compared to other methods, the advantages of our 
method lie in several areas. Firstly, the remarkable proper- 
ties of the shapelet basis functions turn operations such as 
convolution and shear into simple and analytic matrix op- 
erations. In particular, the shear operator can be expressed 
as a simple combination of raising and lowering operators, 
borrowed from the formalism of the QHO. Secondly, it is 
linear in the galaxy intensity, thus avoiding biases intro- 
duced by imperfect knowledge of the noise property of the 
image. Thirdly, the shapelet formalism is capable of using all 
the available shape information, and of optimally estimating 
shear from it. For instance, the KSB method only considers 
gaussian-weighted quadrupole moments, which are exactly 
equal to our second-order shapelet moments. Our method 
thus, in a sense, generalises the KSB approach to include all 
available high-order moments. Finally, the method is ana- 
lytic and well-defined mathematically, and is thus more re- 
liable and stable than KSB which is known to suffer from 
ill-defined quantities (see Kuijken 1999; Kaiser 2000). 



Thanks to these advantages and to its overall complete- 
ness and clarity, our method is well-placed for the analysis of 
current and future weak shear surveys. In particular, it will 
allow us to fully exploit the remarkable resolution of future 
space-based weak lensing surveys with HST and the planned 
SNAP mission (Perlmutter et al. 2001). It thus promises to 
provide the sophistication necessary to enter the next stage 
of high-precision shear analysis. 
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